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We first review how to determine the rate of vibrational energy relaxation (VER) using 
perturbation theory. We then apply those theoretical results to the problem of VER of a 
, CD stretching mode in the protein cytochrome c. We model cytochrome c in vacuum as a 

normal mode system with the lowest-order anharmonic coupling elements. We find that, 
for the "lifetime" width parameter 7 = 3 ~ 30 cm~^, the VER time is 0.2 ~ 0.3 ps, which 
agrees rather well with the previous classical calculation using the quantum correction 
factor method, and is consistent with spectroscopic experiments by Romesberg's group. 
We decompose the VER rate into separate contributions from two modes, and find that 



(N 



PR 

O ! the most significant contribution, which depends on the "lifetime" width parameter, comes 

' ^ ' from those modes most resonant with the CD vibrational mode. 

,cr: 

1 Introduction 

m ' 

I Vibrational energy relaxation (VER) is fundamentally important to our understanding of chem- 

ical reaction dynamics as it influences reaction rates significantly In general, estimat- 
. ing VER rates for selected modes in large molecules is a challenging problem because large 

I molecules involve many degrees of freedom and, furthermore, quantum effects cannot be ig- 

nored [2j. If we assume a weak interaction between the "system" and the surrounding "bath", 
. however, we can derive an estimate of the VER rate through Fermi's golden rule [31l3||SJini: A 

Q I VER rate is written as a Fourier transformation of a force-force correlation function. Though 

it is not trivial to define and justify a separation of a system and a bath, such a formulation 
has been successfully applied to many VER processes in liquids PTI and in proteins ISj- 

Here we apply such theories of VER to the problem of estimating the vibrational population 
relaxation time of a CD stretching mode, in short, a CD mode, in the protein cytochrome c 
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^ ■ [0]. (We will define the CD mode to be the system and the remainder of the protein to be 

the bath.) Recently Romesberg's group succeeded in selectively deuterating a terimnal methyl 
group of a methionie residue in cytochrome c jlOj . The resulting CD mode has a frequency 
ujs — 2100 cm^^, which is located in a transparent region of the density of states of the protein. 
As such, spectroscopic detection of this mode provides clear evidence of the protein dynamics, 
including the VER of the CD vibrational mode. Note that at room temperature (T = 300 K) 
Phujs — 10 where f3 = l/{kBT), hence quantum effects are not negligible for this mode. 
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in mitochondrial inner-membranes, chloroplasts of plants, and bacteria jllj . Its functions are 
related to cell respiration [l^, and cyt c, using its heme molecule, "delivers" an electron from 
cytochrome be 1 to cytochrome oxidase - two larger proteins both embedded in a membrane. 
Recently it was also found that cyt c is released when apoptosis occurs ^Sl- In this sense, cyt 
c governs the "life and death" of a cell. 

The heme molecule in cyt c has a large oscillator strength, and serves as a good optical 
probe. As a result, many spectroscopic experiments have been designed to clarify VER and the 
(un)folding properties of cyt c Cyt c is often employed in numerical simulations jl5l 
because a high resolution structure was obtained ^21 and its simulation has become feasible. 
Attempts have also been made to characterize cyt c through ab initio (DFT) calculations 

[iHiiin]. 

VER of the selected CD mode in the terminal methyl group of methionine (MetSO) was 
previously addressed by Bu and Straub :9 : They used equilibrium simulations for cyt c in 
water with the quantum-correction factor (QCF) method [20], and predicted that the VER 
time for the CD mode is on the order of 0.3 ps. However, their results are approximate: The 
use of the QCF method is not justified a priori, and their analysis is based on a harmonic model 
for cyt c. To extend that previous analysis, in this work, we model cyt c in vacuum as a normal 
mode system and include the lowest anharmonic coupling elements. A similar analysis has 
been completed for another protein myoglobin by Kidera's group |2J and by Leitner's group 
[221 • Use of a reduced model Hamiltonian allows us to investigate the VER rate of the CD 
mode in cyt c more "exactly" and to move beyond the use of quantum correction factors and 
the harmonic approximation. 

This paper is organized as follows: In Sec. [21 we derive the principle VER formula employed 
in our work, and mention the related Maradudin-Fein formula. In Sec. [HI we apply those 
theoretical results for the rate of VER to the CD mode in cyt c, and compare our results with 
the classical simulation by Bu and Straub, and the experiments by Romesberg's group. In 
Sec. ^ we provide a summary of our results, and discuss further aspects of VER processes in 
proteins. 



2 Vibrational energy relaxation (VER) 

2.1 Perturbation expansion for the interaction 

We begin with the von Neumann equation for the complete system written as 

inj^p{t) = [n,p{t)]. (1) 

The interaction representation for the von Neumann equation is 

inj^~p{t) = [v{t)rp{t)i (2) 

where 

W = Wo + V = W^ + Hs + V (3) 

and 

p(t) = e*^"*/^/)(t)e-^^o*/^, V{t) = e^not/hy^-iHot/n^ ^^-^ 

Here TCs is the system Hamiltonian representing a vibrational mode, TCb the bath Hamiltonian 
representing solvent or environmental degrees of freedom, and V the interaction Hamiltonian 
describing the coupling between the system and the bath. An operator with a tilde means the 
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the perturbation expansion for V as 
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m = pm + ^ L dmt'),p{t')] 







Let us calculate the following probability: 

P,{t) ^ TT{p,p{t)} = TT{p,e-''>^^'/''p{ty'^''^''}, (6) 

Pv = \v){v\'^1b, (7) 

p(0) = PS ^PB = \vo) {vol ^e-'^'^'^/ZB, (8) 

Zb = TtbIc-^'^^}, (9) 

where the initial state is assumed to be a direct product state of ps = \vo){vo\ and ps = 
e~^^^ /Zb- Here \v) is the vibrational eigenstate for the system Hamiltonian Tls, i.e., 7{s\v) = 
E^\v). The VER rate Ty^ ^y may be defined as follows: 

r^o^^= lim ^P^(t). (10) 

t^oo at 

Note that the results derived from this definition are equivalent to those derived from Fermi's 
golden rule [22] . Hence we refer to them as a Fermi's golden rule formula. 



2.2 General formula for VER 

First we notice that 

P,{t) = Tr{p,e-^^°*/V(t)e^^«*/^} = Tr{p,p(t)} (11) 

as pv commutes with TLq. If we assume that v / uq, then Pvp{0) = 0. Using this fact, we obtain 
the lowest (second) order result 

Py{t) ^ -L^ fdt' f dt''Tr{p4V{t'),[V{t''),pm]} 
[inY Jo Jo 

= 4 f'dt' r dt''Tr{p,V{t')p{0)V{t'')+p,V{t'')p{0)V{t')} 
n Jo Jo 

^ f dt' f dt"[e''^-o-(t'-nc{t' - t") + e*^"o4i"-t')c(t" - t')\ (12) 

where 



fi^ Jo Jo 



C{t) = (Ko.WV..o(0)) =TrB{pBV.o.WK.o(0)}, (13) 
K,o(t) = {v\V{t)\vo), (14) 
= (£^„o - Ey)/n. (15) 

Hence the lowest order estimate of the VER rate is given by 

r^o^^ = lim 4j /* (it"[e^'^-o-{*-*")c7(t _ t") + e^^"o4«"-i)c(t" _ t)\ 

t^oo fi^ Jo 
1 

= — dte"^^o-tC{t). (16) 

ft J ~^oo 
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expressed as 

V = -qH{<lk],{Pk]) (17) 

where {qk}-,{Pk} are position and momentum variables for the bath. This J^{{qk}-,{Pk}) is a 
force apphed to the system by the bath. Thus we finally obtain the following Fermi's golden 
rule formula El El 

r„„^, = r dte'^^^^\T{t)Tm (18) 

Tl J —CO 

where q^^^ = {vo\q\v), T{t) = (AHBt/nj:^~iHBt/n^ {r{t)T{<d)) = TiB{pB^{t)T{<d)}. In 
most situations, the transition from fo = 1 to w = is considered. In such a case, qio = 
^Jh/2msoJs-, where ms is the system mass and oog = wio is the system frequency in the 
harmonic approximation. Hence 

2.3 Use of a symmetrized auto-correlation function 

It is useful to define a symmetrized force- force correlation function as [HI IS El El 

s{t) = \mt)m) + {mHm- m 

Since S{t) is real and symmetric with respect to t, S{t) = S{—t), we consider it to be analogous 
to Sc\{t), the classical limit of the correlation function. Hereafter we drop the tilde on 
for simplicity. By half-Fourier transforming S{t) with the use of the relation {J^(0)J^{t)) = 
(J^(t - il3h)J='{0)), we have 

/ dte'^^^Sit) = - dte'^^{J^{t)J^{0)) + - dte*^*(.F(0).F(t)) 
Jo 2 Jo 2 Jo 

1 roo 1 roo 

= - dte'^\j^{t)J^iO)) + - dte''^\j^{t-i(3h)J^{0)) 

2 JO 2 Jo 

1 roo 1 rao 

= dte'^\j^it)J^{0)) + ^e-^^^ dte'^\T{t)T{<d)) 

2 Jo 2 Jo 

-{l + e-^^) / dte'''\T{t)T{{))) (21) 
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Taking the real parts of both sides, we have 

/•OO 1 POO 

/ dt cos{tot) S{t) = -il + e-^^^) dte'^\T{t)T{<d)) (22) 
Jo 4 J_oo 

where we have used the fact that S{-t) = S{t) is real and (J^(t)J^(O))* = (J^(O)J^(t)) = 
{T{—t)J^{0)). Hence, Eq. (|19j) can be rewritten as P 

1 2 

ri_o = . / dt cosiujst) S{t). (23) 

msnuJs 1 + e '^"■^s Jq 

Note that this expression diverges in the classical limit because ri_>o oc l/fi. According to 
Bader-Berne 0], to make contact with the classical limit, we introduce another VER rate as 
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±. = (i-e-^^-s)ri^o 



2C(/3, Tiujs) / dt cos{uJst) S{t) (24) 
Jo 
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1 1 _ p-priujs 

C{P, huos) = -r — (25) 

^ ' msfiuJs I + e-f^'^'^s ^ ^ 

This is a final quantum expression, which can be interpreted as an energy relaxation rate and 
be used to estimate the VER rate.^ 

2.4 Quantum correction factor method and other methods 

Though Eq. H24() is exact in a perturbative sense, it is demanding to calculate the quantum 
mechanical force autocorrelation function S{t) even for small molecular systems. Hence, many 
computational schemes have been developed to approximate the quantum mechanical force 
autocorrelation function. 

Skinner and coworkers advocated to use the quantum correction factor (QCF) method pHj . 
which is the replacement of the above formula Eq. ()18() with 

ry^^ = Q{ujs) — / dt cos{ujv„.at) Sci{t) (26) 
n JO 

where Sc\{t) = (J^(t)J^(0))ci and the bracket means a classical ensemble average (not a quantum 
mechanical average). This approach is very intuitive and easily applicable for large molecular 
systems because one only needs to calculate the classical force autocorrelation function Sc\{t) 
multiplied by an appropriate QCF Q{los)- There exist several QCFs corresponding to different 
VER processes |2()] . 

However, one challenge that arises in the application of the QCF method is that we do not 
know a priori which VER process is dominant for the system considered. Furthermore, it is 
possible that several VER processes compete Hence one must be careful in the application 
of the QCF method, and Skinner and coworkers have provided a number of examples of how 
this can be accomplished. 

In this paper, we employ the harmonic approximation for the relaxing oscillator, and the 
vibrational relaxation time 7"^*-'^. Hence Eq. (|26)) is transformed to 

^('^^^ f'"dtcosi.st)Sdt) = ^^, (27) 



^QCF ^^f^^^ V ^ y CIV ; ^^^^ ^ci 

where we have introduced the classical VER rate 

1 B f°° 

T^^^msJQ ^o^^'^^*) '^d(i) (28) 

which is the classical limit 7i ^ of Eq. 1)24(1 . This result can also be derived from a classical 
theory of Brownian motion, and is known as the Landau- Teller-Zwanzig (LTZ) formula. 

Alternatively, one may calculates S{t) itee// systematically using controlled approximations. 
Calculating a correlation function for large systems has a long history in chemical physics |25j . 
including recent applications to VER processes in liquid |26| I27j . The vibrational self-consistent 
field (VSCF) method 28^ will also be useful in this respect. 

On the other hand, if we approximate T as a. simple function of {^fc}, {pfc}, we can calculate 
S{t) rather easily and, in a sense, more "exactly." In the next section, we explore such an 
approach. 

^Although the experimental observable is Fi^o, we note that 1/Ti ~ ri_o because fihus 2> 1 for our case 
of a CD stretching mode. 
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2.5.1 Taylor expansion of the force 

We can formally Taylor-expand the force as a function of the bath variables {qk}, {Pk}' 

HiQk}, {Pk}) = E 4^ Ik + E Bk^Pk + E A^l'lkqk' + E Bj^l'PkPk' + ■■■ (29) 

where the expansion is often truncated in the literature following the first term. Depending on 
the system-bath interaction considered, higher order coupling including the third and fourth 
terms can be relevant. For example, the fourth term appears in benzene to represent the 
interaction between the CH stretch and CCH wagging motion 29 through the Wilson G 
matrix |3nj. In the case of a CD stretching mode in cyt c, as discussed below, or a CN~ 
stretching mode in water .24j, the third term is relevant for VER. 



2.5.2 Contribution from the first term 

If the fir 
becomes 



If the first term J2k'^k^\k is dominant in the force, then the force- force correlation function 



{:F{t)m) = T.4^4^(<ik{t)qk'm = T.^^^ii^^ (30) 
kk' k 

where we have used 



9fcW = y^(afee— ^■* + 4e^"^-*) (31) 
and {a\ak') = nk6k,k' with = l/(e^^'^'= — 1) because 

PB oc e-^«« = e-/5E.'*-fc(«^fc+i/2)_ (32) 
Here we have assumed that the bath Hamiltonian is an ensemble of harmonic oscillators: 

7^B = E ^M4^k + 1/2) = E ( Y + Y*?!) (33) 
where u>k is the k-th mode frequency for the bath, and 



Pk = -i\h^iak-al). (34) 



Thus we obtain 



•5(i) = E o' (2™fe + 1) cos ""kt (35) 
k ^"^k 

and 

^ = TThCiP, hus) E ^^^(2nfc + iMoos -uJk) + S{ujs + u^t)]- (36) 

Tl ^ LOk 

The contribution from the second term Yl,k ^^^Pk is calculated in the same way. 



6 



(2) 

If the third term J2k,k' k'IkQk' is dominant in the force, then 

(mm) = E 4%4'',fc"'(9fe(*)9fc'(i)9fe"(o)9fc"'(o)) 

k,k',k",k"' 

= R—{t) + R++{t) + R+-{t) 



where 



k,k'k''k" 



R+-{t) = ^ E ^S^fe'^^fc''' [^"^4'(4''«fe''' + «fc''4'''))e 
+ (4afc'(4"«fc"' +«ik"4"'))e'^'"="'^'='^* 



with 



.(2) 4(2) 
^(2) _ ^k,k'^k",k"' 

^UJk^k'^k"^k"' 



Using the following 

{akak'a\iia\i,i) = (1 + rafc)(l + nk'){5kk"5k'k"' + ^kW'^k'k 

{a\a\iak"akiii) = nknk'{6kk"^k'k"' + Skk"'Sk'k"), 

{akal,{al„ak"' + ak"al,„)) = {I + nk){l + 2nk")5kk'Sk"k"' 

+ (1 + nk)nk'{6kk"Sk'k"' + ^kk"'^k'k"), 

{alak'{al„ak"' + afe//a^,„)) = nk{l + 2nk")Skk'Sk"k"' 

+nk{l + nk'){Skk"Sk'k"' + Skk"'Sk'k"), 



we have 



^2 



k,k' 

^++\'^) — 2 Z^^k,k',k,k'^knk'e , 
fe,fe' 

fe,fe' 



where we have used Afl = ^g)^, and 



^2 



Hence we obtain 



(•^(0))' = ^ E + 2nfc)(l + 2nfc0. 



^(*) = E [Cfel' cos(a;fc + iOk')t + (^J, cos(a;fc - ujk')t\ + (^(0))^ 



^ k,k' 

+Ck^k'[^(.^k - Wfc' - UJs) + 6{ujk- uJk' + ujs)]\ (51) 
where we have assumed 7^ and defined 



^2 



Ck~k' = -jD'^kl',k,k'i'^k + nk' + 2nknk>). (53) 

Though its appearance is rather different, the formula Eq. (|51|) is equivalent to that derived by 
Kenkre, Tokmakoff, and Fayer |5j as well as by Shiga and Okazaki [2^. There is also a similar 
result known as the Maradudin-Fein formula 

W = W^decay + VFcoU, (54) 

W^decay = TT^^^— "^-^^^(l + J^fc + nfc')'5(ws " ^^fc " t*^fe')> (55) 
2msujs f^, ^k^k' 

Wcow = -^^^y2^^^^^(^k-nk')S{uJs + uJk-uJk'), (56) 
msujs f^, ^k^k' 

which has been utilized to describe VER processes in glasses 32^ and in proteins by Leitner's 
group j22j . As was demonstrated in 5 , this formula is also equivalent to Eq. ()51() : in the 
following we make use of Eq. (|51|). A problem with this formula is that we cannot take its 
continuum limit in the case of finite systems like proteins. As a remedy, a width parameter 
related to the vibrational lifetime is usually introduced leading to a definite value for the VER 
rate. We will discuss this problem in Sec. 13.31 



3 Application to a CD stretching mode in cytochrome c 
3.1 Definition of system and bath 

We take horse heart cytochrome c (cyt c) as an example of how one may estimate the rate of 
VER for selected modes in proteins. We use the CHARMM program j,33j to describe the force 
field, to minimize the structure, and to calculate the normal modes for the system. Starting 
from the IHRC structure for cyt c in Protein Data Bank (PDB) one hydrogen atom of 
the terminal methyl group of MetSO was deuterated. The energy of the protein structure was 
minimized in vacuum using the conjugate gradient algorithm. We diagonalized the hessian 
matrix (second derivatives of the potential) for that mechanically stable configuration of the 
protein: 

dxidxj ^rriimj dxidxj 

where Vcharmm is the CHARMM potential, and Xi = ^Jmlxi are mass weighted Cartesian 
coordinates. The number of atoms in cyt c is 1745 (myoglobin has 2475 atoms), so the hessian 
matrix is 5235 x 5235, and its diagonalization was readily accomplished using the vibran 
facility in CHARMM ^ . 

The result of this calculation was the density of states (DOS) for the system as shown in 
Fig. n The DOS consists of three regions: (1) below around 1700 cm~^, (2) from around 
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to rotational and torsional motions of the protein, and the third to bond stretching motions 
such as CH bonds. The second is rather "transparent" but one can observe one mode localized 
around the CD bond stretching mode in MetSO with frequency 2129.1 cm~^ as shown in Fig. [2 
(a). Hence we refer to this as a CD stretching mode, or CD mode; the dynamics of which is 
the focus of our study. 

0.01 F ^ i ^ q 




1e-005 ' ^ ^ — I 

1000 2000 3000 4000 

0) (cm"^) 



Figure 1: Density of states for cytochrome c in vacuum. 




Figure 2: Normal modes of cytochrome c in vacuum, (a) 4357th mode (CD mode) with 
io = 2129.1 cm-i, (b) 3330th mode with lo = 1330.9 cm-\ (c) 1996th mode with uo = 829.9 
cm^^. Only vectors on the terminal methyl group of MetSO in cyt c are depicted. 

At this level of description, the system is an ensemble of harmonic oscillators, i.e., normal 
modes. Since we are interested in VER of the CD mode, we represent it as a system 

2 2 

n_, _ PCD I ^CD „2 /t;oA 
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= E ( Y + Y'^l) • (59) 

The interaction between the system and bath is described by the interaction Hamiltonian 

V = Wcytc -Hs-Hb (60) 

where TCcytc is the Hamiltonian for the full cyt c protein. We will discuss the content of V in 
the following section. 

3.2 Calculation of the coupling constants 

As in Eq. (|17)1. we assume that the interaction Hamiltonian is of the form 

V = -qcDf" (61) 

and Taylor expand the force as Eq. H29() . The first and second terms do not appear because this 
is a normal mode expansion, and the fourth term does not appear as the original coordinates 
are Cartesian coordinates. As in the first approximation, we take the force to be 



•^ = E^fcl'9fc9fc'. (62) 



^(2) 



(2) 

The coupling coefficients A), j are calculated as 

a(^) - ^ (63) 

A problem arises: How does one calculate these coupling coefficients? The most direct approach 
is to use a finite difference method: 

A2) ^ 1 - V-++ - - + V~~+ + V~+- + V+- - V^- 

''^ 2 (2AgcD)(2A(?fc)(2A(?0 ^ > 

where V±±± = y(ibAg'c'£), ztAg^, ibAg;). However, this is rather cumbersome. Instead, we use 
the approximation |221 124j : 

(65) 

where Uik is an orthogonal matrix that diagonalizes the hessian matrix at the mechanically 
stable structure Kij, and Kij{-^IS.qcD) is a hessian matrix calculated at a shifted structure 
along the direction of the CD mode with a shift ^lS.qcD- This expression is approximate but 
readily implemented using the CHARMM facility to compute the hessian matrix. A comparison 
between Eqs. H64|) and H65|) is made in Tabled We also examined the convergence of the results 
by changing /S.qcD-, and found that Aqcn = 0.02A is sufficient for the following calculations. 

The numerical results for the coupling elements are shown in Fig. 13 The histogram for the 
elements is shown in Fig. |1J As one can see from these figures, most of the elements are small, 
while the largest coupling elements are rather large. See Table El Note that the combination 
(3330, 1996) is particularly significant for the CD mode because it approximately satisfies the 
resonant condition (21] : 

|u;cD-c^fc-u;,| «0(|4'J|). (66) 

As shown in Fig.|51(b), (c), these modes are localized near the terminal methyl group of MetSO 
as well as the CD mode. In such a case, resonant energy transfer (Fermi resonance) is expected 
as shown by Moritsugu-Miyashita-Kidera |2J|- We have observed similar behavior in cyt c 
when the CD mode was excited and the energy immigration to other normal modes facilitated 
by resonance was followed. 
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Table 1: Comparison between the finite difference method Eq. H64() and the formula Eq. ()65() . 
We have used AqcD = 0.02 A, and A^^j is given in kcal/mol/A'^. Note that {k,l) are mode 
numbers, not wavenumbers. 



{k,l) 


formula, Eq. (j65j) 


Finite difference 


(3330, 1996) 


22.3 


22.4 


(3330, 4399) 


-29.6 


-29.5 


(3327, 1996) 


-5.7 


-5.8 


(1996, 678) 


0.64 


0.63 



Table 2: The largest coupling elements. The value of A^^j are given in kcal/mol/A^. Note 
that (k, I) are mode numbers, not wavenumbers. 



ik,l) 


\^k,i 1 


(1996, 1996) 


42.9 


(4399, 3330) 


29.6 


(4622, 3170) 


27.3 


(3330, 1996) 


22.3 
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k 



1000 2000 3000 4000 5000 
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(2) 

Figure 3: Distribution of the coupling elements. Left: 50.0 > |^^/| > 5.0, Right: 5.0 > 

\A^^j\ > 0.5. The value of A^^j are given in units of kcal/mol/A^. Note that {k,l) are mode 
numbers, not wavenumbers. 



11 




IASI (kcal/mol/A ) 



Figure 4: Histogram for the amplitude of the couphng elements. 
3.3 Assignment of the "lifetime" parameter 

We cannot directly evaluate Eq. (|51() as it contains delta functions. Evaluation of this expression 
for a finite system like a protein leads to a null result. To circumvent this problem, we "thaw" 
the delta function 5{x) as 



6{x) = -- 
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2^ 2 (67) 

using a width parameter 7. Physically this means that each normal mode should have a lifetime 
~ 1/7 due to coupling to other degrees of freedom, i.e., the surrounding environment including 
water (or we might be able to interpret I/7 as a time resolution) . It is difficult to derive 7 from 
first principles, so we take it to be a phenomenological parameter as in the literature [221 132j . 
As a result, the VER rate, Eq. (|5T|) . for the CD mode becomes 



^ k,k' 



+ 



+ 



l'^k,k' 



l^k,k' 



+ 



^CkJ' 



72 + {uk - iOk' - ^^s)'^ 7^ + i^k - ^k> + ^s) 



(68) 



We employ this expression in our subsequent calculations. 



3.4 Results 

3.4.1 Classical calculation 

Classical force-force correlation functions and their (cosine) Fourier transformations are shown 
in the left and middle column of Fig. [3 for five different trajectories. Here we have defined a C, 
function as 

m = -^S,iit), (69) 
ms 
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Figure 5: Left: Classical data for the force-force correlation function. Middle: Fourier spectra 
for the correlation function. Right: The corresponding coarse-grained Fourier spectra. The 
"lifetime" width parameter 7 = 3 cm"-*-. 



13 



are obtained from molecular dynamics simulations of cyt c in water 0. As can be seen, the 
correlation functions oscillate wildly, and the (cosine) Fourier transformations are messy. As 
such, it is difficult to extract a reliable and stable value for the VER rate. 
To address this problem, we introduce the window function 

= exp(-7t). (70) 

The C functions are multiplied by this function and (cosine) Fourier transformed. This corre- 
sponds to broadening each peak of a spectrum with a Lorentzian with width 7. The results 
for five trajectories are shown in the right column of Fig. [3 (The width parameter is taken 
as 7 = 3 cm~^.) The results in the right column are better behaved than those in the middle 
column but there still remain some fluctuations. 

According to Bu and Straub simulations of cyt c in water I^J, we take cos = 2135 cm~^ 
to investigate the 7 dependence of the result as shown in the left of Fig. El We see that 
C{u!s) — 1-1 ~ 1-2 ps~^ for 7 ~ 3 ~ 30 cm~^. Since Q{ujs)/{P^ujs) — 2.4 ~ 3.0 for two-phonon 
processes jH], this corresponds to a VER time of 0.3 ~ 0.4 ps according to Eq. ((^.^ 



3.4.2 Quantum calculation 

We use the formula Eq. ()68|) as a quantum mechanical estimate of the VER rate. The 7 
dependence of the result is shown on the right hand side of Fig. El We see that, for 7 ~ 3 ~ 30 
cm~^, the quantum mechanical estimate gives Ti ~ 0.2 ~ 0.3 ps, which is similar to the classical 
estimate Eq. Ti ~ 0.3 ~ 0.4 ps. 




0.01 



0.1 



1 10 
7(cm ') 



100 



1000 



0.01 



0.1 



1 10 
7(cm ') 



100 1000 



Figure 6: Left: Classical VER rate for five trajectories as a function of the "lifetime" width 
parameter 7. Right: VER rate calculated by Eq. (|68]) as a function of 7. 

In Tables O we list the largest contributions to the VER rate for different width parameters. 
For the case of 7 = 3 cm~^, the largest contribution is due to modes (3823,1655). This 
combination of modes is nearly resonant with the CD mode as |<-i;3823+i^1655~'-^cd| ~ 0.03 cm~^. 
Though the coupling element for the combination is small (1^3323 1655I =5.1 kcal/mol/A^), this 
mode combination contributes significantly to the VER rate. On the other hand, for the case of 
7 = 30 cm~^, the largest contribution results from the combination of modes (3330,1996). This 
combination is somewhat off-resonant, i.e., |i-i;333o + wiqqq — locdI = 32 cm~^, but the coupling 
element is very large (1^3330 iggel = 22.3 kcal/mol/A^), and the contribution is significant. 

^In the VER calculation of myoglobin |22|. Leitner's group took 7 = 0.5 ~ 10 cm^^ to be the width, and 
confirmed that the result is relatively insensitive to the choice of 7 in this range. 
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non-negligible contributions from other combinations of modes. 

Table 3: The largest contributions to the VER rate (in units of ps~^) for 7 = 3 cm~^ (left) 
and 7 = 30 cm~^ (right). Note that {k,l) are mode numbers, not wavenumbers. 



{k, I) 


contribution 


{k,l) 


contribution 


(3823, 1655) 


1.10 ( 19 %) 


(3330, 1996) 


0.88 (22 %) 


(3823, 1654) 


0.43 ( 8 %) 


(3823, 1655) 


0.11 (3 %) 


(3822, 1655) 


0.37 ( 6 %) 


(3170, 2196) 


0.07 (2 %) 


(3330, 1996) 


0.17 ( 3 %) 


(1996, 1996) 


0.05 (1 %) 


(3822, 1654) 


0.15 ( 3 %) 


(3823, 1654) 


0.04 (1 %) 


(3823, 1661) 


0.14 ( 3 %) 


(3170, 2202) 


0.04 (1 %) 


(3822, 1661) 


0.05 ( 1 %) 


(3822, 1655) 


0.04 (1 %) 


(3822, 1656) 


0.05 ( 1 %) 


(3327, 1996) 


0.03 (1 %) 


(3823, 1658) 


0.04 ( 1 %) 


(3330, 1655) 


0.02 (1 %) 



We have also examined the temperature dependence of the VER rates using Eq. (|68|) . As 
shown in Fig. d for T < 300 K there is little temperature dependence as has been addressed 
in the case of myoglobin |^. Thus we can say that the relaxation of the CD mode is quantum 
mechanical rather than thermal because the decay at 300 K is similar to that at 0.3 K. 
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Figure 7: Temperature dependence of the VER rate calculated by Eq. (|68|) . The width 
parameter is 7 = 3 cm~^. 

3.4.3 Discussion 

We examine the relationship between the theoretical results described above and the corre- 
sponding experiments of Romesberg's group which has studied the spectroscopic properties of 
the CD mode in cyt c TD . They measured the shifts and widths of the spectra for different 
forms of cyt c; the widths of the spectra (FWHM) were found to be Awfwhm — 6.0 ~ 13.0 
cm~^. A rough estimate of the VER rate leads to 

Ti ~ 5.3/A6JFWHM (ps) (71) 
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prediction computed using Eq. and appropriate QCFs (0.3 ~ 0.4 ps) and the perturbative 
quantum mechanical estimate using the reduced model (0.2 ~ 0.3 ps). This result appears to 
justify the use of QCFs and the reduced model in this situation, and suggests that the effects 
of the protein solvation (by water) are negligible in describing the VER of the CD mode. 
Of course, we must be careful in comparing the estimate derived from (jTlf) as there may be 
inhomogeneity in the experimental spectra.^ As such, it is more desirable to calculate not only 
VER rates but spectroscopic observables themselves to compare with experiments. 

Finally we discuss the relation between this work and previous work on carbon-monoxide 
myoglobin (MbCO). Though there are many experimental studies on Mb jHEl) we focus on the 
experiments of Anfinrud's group [3^] and Fayer's group |37| on MbCO. The former group found 
that the VER time for CO in the heme pocket (photolyzed MbCO) is ~ 600 ps, whereas the 
latter group found that the VER time for CO bounded to the heme is ~ 20 ps. This difference 
is interpreted as follows: CO is covalently bonded to the heme for the latter case, whereas CO 
is "floating" in the pocket in the former case, i.e., the force applied to CO for the latter case is 
stronger than that for the former case. This difference in the magnitude of the force causes the 
slower VER for the "floating" CO. In this respect, the CD bond is expected to be stronger than 
the CO-heme coupling. This may explain a VER time ~ 0.1 ps, which is similar to the VER 
times for the CH(CD) stretching modes in benzene (or perdeuterobenzene) |29| I38j . It will be 
interesting to apply a similar reduced model to the analysis of the VER of CO in MbCO. 

4 Summary and further aspects 

After reviewing the VER rate formula derived from quantum mechanical perturbation calcu- 
lations, we applied it to the analysis of VER of a CD stretching mode in cyt c. We modeled 
cyt c in vacuum as a normal mode system with the third-order anharmonic coupling elements, 
which were calculated from the CHARMM potential. We found that, for the width parameter 
7 = 3 30 cm~^, the VER time is 0.2 ~ 0.3 ps, which agrees rather well with the previous 
classical calculation using the quantum correction factor (QCF) method, and is consistent with 
the experiments by Romesberg's group. This result indicates that the use of QCFs or a reduced 
model Hamiltonian can be justified a posteriori to describe the VER problem. We decomposed 
the VER rate into contributions from two modes, and found that the most significant contri- 
bution, which depends on the "lifetime" width parameter, results from modes most resonant 
with the CD mode. 

Finally we note several future directions which should be studied: (a) Our final results for 
the VER rate depend on a width parameter 7. Unfortunately we do not know which value 
is the most appropriate for 7. Non-equilibrium simulations (with some quantum corrections 
pUj ) might help this situation, and are useful to investigate energy path ways or sequential 
IVR (intramolecular vibrational energy redistribution) jlHj in a protein, (b) This work is 
motivated by pioneering spectroscopic experiments by Romesberg's group. The calculation of 
the VER rate and the linear or nonlinear response functions, related to absorption or 2D-IR 
(or 2D-Raman) spectra ESI 113) is desirable, (c) Romesberg's group investigated a 

spectroscopic change due to the oxidation or reduction of Fe in the heme; such an electron 
transfer process t45_ is fundamental for the functionality of cyt c. To survey this process 
dynamically, it will be necessary to combine some quantum chemistry (ab initio) calculations 
with MD simulations [TSl 11511^ . 

The authors thank Dr. J. Gong for noting Ref. 0, Prof. A. Kidera, Prof. S. Okazaki, 
Prof. J. L. Skinner, Prof. D.M. Leitner, Prof. Y. Mizutani, Dr. T. Takami, Dr. T. Miyadera, 

■^We have confirmed that the methyl group does not rotate during the equihbrium simulations. Thus we can 
exclude the rotation as a possible reason of inhomogeneity. 
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